Fit models to biomass density and predict covariates
Max Lindmark
07 November, 2022
TODO
- Fix pred grid and standardization to account for different sizes!
Intro
In this script, I take the collated stomach data set and calculate aggregates (feeding ratio, total weight of prey groups) and predictor variables for diet data, aggregate to get 1 stomach = 1 row per prey type (not prey individual). I also select only the columns I need for model fitting, join environmental covariates and cpue covariates for cod and flounder, and lastly saduria biomass densities.
Load packages & source functions
# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(readxl)
library(tidylog)
library(RCurl)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(patchwork)
library(janitor)
library(forcats)
library(gapminder)
library(viridis)
library(ggridges)
library(raster)
library(icesDatras)
library(ggalluvial)
library(ggrepel)
library(ncdf4)
library(chron)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(geosphere)
library(quantreg)
#> Warning in .recacheSubclasses(def@className, def, env): undefined subclass
#> "numericVector" of class "Mnumeric"; definition not updated
library(brms)
#> Warning: package 'Rcpp' was built under R version 4.0.5
library(sdmTMB)
options(mc.cores = parallel::detectCores())
world <- ne_countries(scale = "medium", returnclass = "sf")
# Source code for map plots
# Source code for map plots
source("/Users/maxlindmark/Dropbox/Max work/R/spatial_metabolic_index/R/functions/map_plot.R")
plot_map_fc <- plot_map_fc + theme(legend.position = "bottom")
theme_set(theme_plot())
# Continuous colors
options(ggplot2.continuous.colour = "viridis")
# Discrete colors
scale_colour_discrete <- function(...) {
scale_colour_brewer(palette = "Set2")
}
scale_fill_discrete <- function(...) {
scale_fill_brewer(palette = "Set2")
}Read data
# Load cache
# qwraps2::lazyload_cache_dir(path = "R/01_sdm_mi_cache/html")
pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
#> Rows: 129346 Columns: 29
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): ices_rect
#> dbl (28): X, Y, depth, year, oxy, temp, lon, lat, sub_div, density_saduria, ...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")
#> Rows: 120107 Columns: 29
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): ices_rect
#> dbl (28): X, Y, depth, year, oxy, temp, lon, lat, sub_div, density_saduria, ...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid <- bind_rows(pred_grid1, pred_grid2)
pred_grid <- pred_grid %>%
dplyr::select(-density_saduria, -biomass_spr, -biomass_her, -biomass_spr_sd, -biomass_her_sd,
-density_cod, -density_fle, -depth_rec, -temp_rec, -oxy_rec, -density_cod_rec,
-density_cod_sd, -density_fle_sd, -density_fle_rec, -density_saduria_rec,
-density_saduria_sd, -temp_sd, -oxy_sd) %>%
drop_na(oxy, temp)
#> drop_na: no rows removed
cpue <- readr::read_csv("data/full_cpue_length_q1_q4.csv") %>%
janitor::clean_names() %>%
rename(X = x,
Y = y) %>%
filter(quarter %in% c(1, 4))
#> Rows: 1759609 Columns: 18
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (6): haul.id, IDx, ices_rect, id_haul_stomach, species, haul.id.size
#> dbl (12): density, year, lat, lon, quarter, sub_div, length_cm, depth, temp,...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> rename: renamed 2 variables (X, Y)
#>
#> filter: removed 21,168 rows (1%), 1,738,441 rows remaining
# Match pred grid else we cannot predict later
cpue <- cpue %>% filter(year %in% unique(pred_grid$year))
#> filter: removed 43,092 rows (2%), 1,695,349 rows remaining
pred_grid
#> # A tibble: 249,453 × 11
#> X Y depth year oxy temp lon lat sub_div ices_rect depth_sd
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
#> 1 306 6012 13.2 1993 4.11 9.78 12.0 54.2 24 37G2 21
#> 2 306 6012 13.2 1994 6.27 9.66 12.0 54.2 24 37G2 21
#> 3 306 6012 13.2 1995 5.42 9.63 12.0 54.2 24 37G2 21
#> 4 306 6012 13.2 1996 6.11 9.99 12.0 54.2 24 37G2 21
#> 5 306 6012 13.2 1997 5.78 9.84 12.0 54.2 24 37G2 21
#> 6 306 6012 13.2 1998 6.45 8.56 12.0 54.2 24 37G2 21
#> 7 306 6012 13.2 1999 5.53 10.2 12.0 54.2 24 37G2 21
#> 8 306 6012 13.2 2000 4.37 12.5 12.0 54.2 24 37G2 21
#> 9 306 6012 13.2 2001 5.81 9.91 12.0 54.2 24 37G2 21
#> 10 306 6012 13.2 2002 5.26 9.92 12.0 54.2 24 37G2 21
#> # … with 249,443 more rows
# Now scale variables in the cpue data by the mean and sd in the prediction grid
# pred_grid$log_depth <- log(pred_grid$depth)
# pred_grid$log_depth2 <- pred_grid$log_depth * pred_grid$log_depth
# mean_depth <- mean(pred_grid$log_depth)
# sd_depth <- sd(pred_grid$log_depth)
# mean_depth2 <- mean(pred_grid$log_depth2)
# sd_depth2 <- sd(pred_grid$log_depth2)
# Calculating a metabolic index
# A_0 is the ratio of constant terms in the Arrhenius equations describing the rates of supply and demand
# B is individual body size (g)
# n is the allometric exponent (difference delta - epsilon)
# E0 is the difference in temperature dependences of supply and demand
# T is temperature in degrees Kelvin
# Tref is an arbitrarily chosen reference temperature (15°C)
# kb is Boltzmann’s constant (eV K−1).
# Parameters from Deutsch et al 2015 Science Supp Fig. S2
A_conc = 9.7*10^-15
E_conc <- 0.72
n <- -0.21
Tref <- 273.15 + 10
kb <- 0.000086173324
pred_grid$T_K <- pred_grid$temp + 273.15
pred_grid$oxy2 <- (pred_grid$oxy * 10^3) / 22.391
pred_grid$phi <- A_conc*(2000^n)*(pred_grid$oxy2 / exp(-E_conc / (kb * pred_grid$T_K)))
# Scale variables in data w.r.t. prediction grid
mean_phi<- mean(pred_grid$phi)
sd_phi <- sd(pred_grid$phi)
mean_depth <- mean(pred_grid$depth)
sd_depth <- sd(pred_grid$depth)
mean_oxy <- mean(pred_grid$oxy)
sd_oxy <- sd(pred_grid$oxy)
mean_temp <- mean(pred_grid$temp)
sd_temp <- sd(pred_grid$temp)
# Calculate phi in data
# Metabolic index in the data
B_small_cod <- 0.01*median(filter(cpue, length_cm <= 30 & species == "cod")$length_cm)^3
#> filter: removed 1,415,326 rows (83%), 280,023 rows remaining
B_large_cod <- 0.01*median(filter(cpue, length_cm > 30 & species == "cod")$length_cm)^3
#> filter: removed 819,148 rows (48%), 876,201 rows remaining
# https://github.com/fate-spatialindicators/SDM_O2/blob/master/code/mi_functions.R
cpue <- cpue %>%
mutate(T_K = temp + 273.15,
oxy2 = oxy * (10^3) / 22.391,
size_cl = ifelse(length_cm > 30, "large", "small"),
phi = ifelse(size_cl == "large",
A_conc*(5000^n)*(oxy2 / exp(-E_conc / (kb * T_K))),
A_conc*(30^n)*(oxy2 / exp(-E_conc / (kb * T_K)))))
#> mutate: new variable 'T_K' (double) with 8,652 unique values and <1% NA
#> new variable 'oxy2' (double) with 8,651 unique values and <1% NA
#> new variable 'size_cl' (character) with 2 unique values and 0% NA
#> new variable 'phi' (double) with 17,303 unique values and <1% NA
ggplot(cpue, aes(phi)) + geom_histogram() + facet_wrap(~size_cl, ncol = 1)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 4487 rows containing non-finite values (stat_bin).
cpue <- cpue %>%
mutate(
# log_depth = log(depth),
# log_depth2 = log_depth*log_depth,
# depth_sc = (log_depth - mean_depth) / sd_depth,
# depth_sc_sq = (log_depth2 - mean_depth2) / sd_depth2,
depth_sc = (depth - mean_depth) / sd_depth,
temp_sc = (temp - mean_temp) / sd_temp,
oxy_sc = (oxy - mean_oxy) / sd_oxy,
phi_sc = (phi - mean_phi) / sd_phi)
#> mutate: new variable 'depth_sc' (double) with 150 unique values and 0% NA
#> new variable 'temp_sc' (double) with 8,652 unique values and <1% NA
#> new variable 'oxy_sc' (double) with 8,651 unique values and <1% NA
#> new variable 'phi_sc' (double) with 17,303 unique values and <1% NASummarise data
# Summaries density by length group per haul
cod_cpue <- cpue %>%
filter(species == "cod") %>%
mutate(size_group = ifelse(length_cm >= 30, "large", "small")) %>%
group_by(haul_id, size_group) %>%
summarise(density = sum(density)) %>%
ungroup()
#> filter: removed 539,125 rows (32%), 1,156,224 rows remaining
#> mutate: new variable 'size_group' (character) with 2 unique values and 0% NA
#> group_by: 2 grouping variables (haul_id, size_group)
#> summarise: now 18,066 rows and 3 columns, one group variable remaining (haul_id)
#> ungroup: no grouping variables
cod_cpue_small <- cod_cpue %>%
filter(size_group == "small") %>%
dplyr::select(-size_group) %>%
rename(density_cod_small = density)
#> filter: removed 9,033 rows (50%), 9,033 rows remaining
#> rename: renamed one variable (density_cod_small)
cod_cpue_large <- cod_cpue %>%
filter(size_group == "large") %>%
dplyr::select(-size_group) %>%
rename(density_cod_large = density)
#> filter: removed 9,033 rows (50%), 9,033 rows remaining
#> rename: renamed one variable (density_cod_large)
fle_cpue <- cpue %>%
filter(species == "flounder") %>%
mutate(size_group = ifelse(length_cm >= 22, "large", "small")) %>%
group_by(haul_id, size_group) %>%
summarise(density = sum(density)) %>%
ungroup()
#> filter: removed 1,156,224 rows (68%), 539,125 rows remaining
#> mutate: new variable 'size_group' (character) with 2 unique values and 0% NA
#> group_by: 2 grouping variables (haul_id, size_group)
#> summarise: now 17,622 rows and 3 columns, one group variable remaining (haul_id)
#> ungroup: no grouping variables
fle_cpue_small <- fle_cpue %>%
filter(size_group == "small") %>%
dplyr::select(-size_group) %>%
rename(density_fle_small = density)
#> filter: removed 8,811 rows (50%), 8,811 rows remaining
#> rename: renamed one variable (density_fle_small)
fle_cpue_large <- fle_cpue %>%
filter(size_group == "large") %>%
dplyr::select(-size_group) %>%
rename(density_fle_large = density)
#> filter: removed 8,811 rows (50%), 8,811 rows remaining
#> rename: renamed one variable (density_fle_large)
# Check 0 catches
cod_cpue_small %>%
group_by(haul_id) %>%
summarise(haul_dens = sum(density_cod_small)) %>%
ungroup() %>%
filter(!haul_dens == 0)
#> group_by: one grouping variable (haul_id)
#> summarise: now 9,033 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> filter: removed 1,685 rows (19%), 7,348 rows remaining
#> # A tibble: 7,348 × 2
#> haul_id haul_dens
#> <chr> <dbl>
#> 1 1993:1:GFR:SOL:H20:21:1 15.1
#> 2 1993:1:GFR:SOL:H20:23:31 0.00919
#> 3 1993:1:GFR:SOL:H20:25:2 278.
#> 4 1993:1:GFR:SOL:H20:26:3 112.
#> 5 1993:1:GFR:SOL:H20:27:27 3.64
#> 6 1993:1:GFR:SOL:H20:28:24 741.
#> 7 1993:1:GFR:SOL:H20:29:29 0.174
#> 8 1993:1:GFR:SOL:H20:31:25 20.3
#> 9 1993:1:GFR:SOL:H20:32:5 571.
#> 10 1993:1:GFR:SOL:H20:33:6 282.
#> # … with 7,338 more rows
cod_cpue_large %>%
group_by(haul_id) %>%
summarise(haul_dens = sum(density_cod_large)) %>%
ungroup() %>%
filter(!haul_dens == 0)
#> group_by: one grouping variable (haul_id)
#> summarise: now 9,033 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> filter: removed 1,200 rows (13%), 7,833 rows remaining
#> # A tibble: 7,833 × 2
#> haul_id haul_dens
#> <chr> <dbl>
#> 1 1993:1:GFR:SOL:H20:21:1 142.
#> 2 1993:1:GFR:SOL:H20:22:32 10.3
#> 3 1993:1:GFR:SOL:H20:23:31 5.22
#> 4 1993:1:GFR:SOL:H20:24:30 19.4
#> 5 1993:1:GFR:SOL:H20:25:2 531.
#> 6 1993:1:GFR:SOL:H20:26:3 251.
#> 7 1993:1:GFR:SOL:H20:27:27 19.5
#> 8 1993:1:GFR:SOL:H20:28:24 1561.
#> 9 1993:1:GFR:SOL:H20:29:29 30.6
#> 10 1993:1:GFR:SOL:H20:30:28 6.84
#> # … with 7,823 more rows
fle_cpue_small %>%
group_by(haul_id) %>%
summarise(haul_dens = sum(density_fle_small)) %>%
ungroup() %>%
filter(!haul_dens == 0)
#> group_by: one grouping variable (haul_id)
#> summarise: now 8,811 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> filter: removed 3,236 rows (37%), 5,575 rows remaining
#> # A tibble: 5,575 × 2
#> haul_id haul_dens
#> <chr> <dbl>
#> 1 1993:1:GFR:SOL:H20:21:1 0.554
#> 2 1993:1:GFR:SOL:H20:23:31 0.382
#> 3 1993:1:GFR:SOL:H20:26:3 1.16
#> 4 1993:1:GFR:SOL:H20:27:27 3.50
#> 5 1993:1:GFR:SOL:H20:29:29 0.498
#> 6 1993:1:GFR:SOL:H20:30:28 0.297
#> 7 1993:1:GFR:SOL:H20:31:25 0.526
#> 8 1993:1:GFR:SOL:H20:32:5 0.779
#> 9 1993:1:GFR:SOL:H20:33:6 1.51
#> 10 1993:1:GFR:SOL:H20:34:7 10.1
#> # … with 5,565 more rows
fle_cpue_large %>%
group_by(haul_id) %>%
summarise(haul_dens = sum(density_fle_large)) %>%
ungroup() %>%
filter(!haul_dens == 0)
#> group_by: one grouping variable (haul_id)
#> summarise: now 8,811 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> filter: removed 1,024 rows (12%), 7,787 rows remaining
#> # A tibble: 7,787 × 2
#> haul_id haul_dens
#> <chr> <dbl>
#> 1 1993:1:GFR:SOL:H20:21:1 18.3
#> 2 1993:1:GFR:SOL:H20:22:32 13.3
#> 3 1993:1:GFR:SOL:H20:23:31 1.52
#> 4 1993:1:GFR:SOL:H20:24:30 3.78
#> 5 1993:1:GFR:SOL:H20:25:2 33.4
#> 6 1993:1:GFR:SOL:H20:26:3 25.8
#> 7 1993:1:GFR:SOL:H20:27:27 20.9
#> 8 1993:1:GFR:SOL:H20:28:24 21.2
#> 9 1993:1:GFR:SOL:H20:29:29 7.88
#> 10 1993:1:GFR:SOL:H20:30:28 9.28
#> # … with 7,777 more rowscod_cpue_small <- left_join(cod_cpue_small,
cpue %>%
filter(species == "cod") %>%
distinct(haul_id, .keep_all = TRUE) %>%
dplyr::select(haul_id, year, quarter, X, Y, oxy, temp, depth,
oxy_sc, temp_sc, depth_sc, phi, phi_sc)) %>%
drop_na(phi)
#> filter: removed 539,125 rows (32%), 1,156,224 rows remaining
#> distinct: removed 1,147,191 rows (99%), 9,033 rows remaining
#> Joining, by = "haul_id"
#> left_join: added 12 columns (year, quarter, X, Y, oxy, …)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 9,033
#> > =======
#> > rows total 9,033
#> drop_na: removed 26 rows (<1%), 9,007 rows remaining
cod_cpue_large <- left_join(cod_cpue_large,
cpue %>%
filter(species == "cod") %>%
distinct(haul_id, .keep_all = TRUE) %>%
dplyr::select(haul_id, year, quarter, X, Y, oxy, temp, depth,
oxy_sc, temp_sc, depth_sc, phi, phi_sc)) %>%
drop_na(phi)
#> filter: removed 539,125 rows (32%), 1,156,224 rows remaining
#> distinct: removed 1,147,191 rows (99%), 9,033 rows remaining
#> Joining, by = "haul_id"left_join: added 12 columns (year, quarter, X, Y, oxy, …)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 9,033
#> > =======
#> > rows total 9,033
#> drop_na: removed 26 rows (<1%), 9,007 rows remaining
fle_cpue_small <- left_join(fle_cpue_small,
cpue %>%
filter(species == "flounder") %>%
distinct(haul_id, .keep_all = TRUE) %>%
dplyr::select(haul_id, year, quarter, X, Y, oxy, temp, depth,
oxy_sc, temp_sc, depth_sc))
#> filter: removed 1,156,224 rows (68%), 539,125 rows remaining
#> distinct: removed 530,314 rows (98%), 8,811 rows remaining
#> Joining, by = "haul_id"left_join: added 10 columns (year, quarter, X, Y, oxy, …)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 8,811
#> > =======
#> > rows total 8,811
fle_cpue_large <- left_join(fle_cpue_large,
cpue %>%
filter(species == "flounder") %>%
distinct(haul_id, .keep_all = TRUE) %>%
dplyr::select(haul_id, year, quarter, X, Y, oxy, temp, depth,
oxy_sc, temp_sc, depth_sc))
#> filter: removed 1,156,224 rows (68%), 539,125 rows remaining
#> distinct: removed 530,314 rows (98%), 8,811 rows remaining
#> Joining, by = "haul_id"left_join: added 10 columns (year, quarter, X, Y, oxy, …)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 8,811
#> > =======
#> > rows total 8,811Fit models
# Small cod
spde_small_cod <- make_mesh(cod_cpue_small, xy_cols = c("X", "Y"),
n_knots = 150,
type = "kmeans", seed = 42)
# Large cod
spde_large_cod <- make_mesh(cod_cpue_large, xy_cols = c("X", "Y"),
n_knots = 150,
type = "kmeans", seed = 42)
# Small flounder
spde_small_fle <- make_mesh(fle_cpue_small, xy_cols = c("X", "Y"),
n_knots = 150,
type = "kmeans", seed = 42)
# Large flounder
spde_large_fle <- make_mesh(fle_cpue_large, xy_cols = c("X", "Y"),
n_knots = 150,
type = "kmeans", seed = 42)
plot(spde_small_cod)
plot(spde_large_cod)
plot(spde_small_fle)
plot(spde_large_fle)
Small cod
# Small cod model
mcod_s <- sdmTMB(density_cod_small ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + oxy_sc + temp_sc,
data = cod_cpue_small,
mesh = spde_small_cod,
family = tweedie(link = "log"),
spatiotemporal = "off",
spatial = "on",
time = "year",
reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.056067183534736.
mcod_s_v2 <- run_extra_optimization(mcod_s, nlminb_loops = 1, newton_loops = 1)sanity(mcod_s_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large# Small cod model
mcod_s1 <- sdmTMB(density_cod_small ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + breakpt(oxy_sc) + temp_sc,
data = cod_cpue_small,
mesh = spde_small_cod,
family = tweedie(link = "log"),
spatiotemporal = "off",
spatial = "on",
time = "year",
reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.011874176036468.
mcod_s1_v2 <- run_extra_optimization(mcod_s1, nlminb_loops = 1, newton_loops = 1)sanity(mcod_s1_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large# Small cod model
mcod_s2 <- sdmTMB(density_cod_small ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + oxy_sc * temp_sc,
data = cod_cpue_small,
mesh = spde_small_cod,
family = tweedie(link = "log"),
spatiotemporal = "off",
spatial = "on",
time = "year",
reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0247736933054483.
mcod_s2_v2 <- run_extra_optimization(mcod_s2, nlminb_loops = 1, newton_loops = 1)sanity(mcod_s2_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large# Small cod model
mcod_s3 <- sdmTMB(density_cod_small ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + breakpt(phi_sc),
data = cod_cpue_small,
mesh = spde_small_cod,
family = tweedie(link = "log"),
spatiotemporal = "off",
spatial = "on",
time = "year",
reml = FALSE)
mcod_s3_v2 <- run_extra_optimization(mcod_s3, nlminb_loops = 1, newton_loops = 1)sanity(mcod_s3_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large# Small cod model
mcod_s4 <- sdmTMB(density_cod_small ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + s(phi_sc),
data = cod_cpue_small,
mesh = spde_small_cod,
family = tweedie(link = "log"),
spatiotemporal = "off",
spatial = "on",
time = "year",
reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0109864218642901.
mcod_s4_v2 <- run_extra_optimization(mcod_s4, nlminb_loops = 1, newton_loops = 1)sanity(mcod_s4_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large#mcmc_res_sc <- residuals(mcod_s_v2, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
#qqnorm(mcmc_res_sc);qqline(mcmc_res_sc)Large cod
# Large cod model
mcod_l <- sdmTMB(density_cod_large ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + oxy_sc + temp_sc,
data = cod_cpue_large,
mesh = spde_large_cod,
family = tweedie(link = "log"),
spatiotemporal = "off",
spatial = "on",
time = "year",
reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0333087100107043.
mcod_l_v2 <- run_extra_optimization(mcod_l, nlminb_loops = 1, newton_loops = 1)sanity(mcod_l_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large# Large cod model
mcod_l1 <- sdmTMB(density_cod_large ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + breakpt(oxy_sc) + temp_sc,
data = cod_cpue_large,
mesh = spde_large_cod,
family = tweedie(link = "log"),
spatiotemporal = "off",
spatial = "on",
time = "year",
reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0103058445924128.
mcod_l1_v2 <- run_extra_optimization(mcod_l1, nlminb_loops = 1, newton_loops = 1)sanity(mcod_l1_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large# Large cod model
mcod_l2 <- sdmTMB(density_cod_large ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + oxy_sc*temp_sc,
data = cod_cpue_large,
mesh = spde_large_cod,
family = tweedie(link = "log"),
spatiotemporal = "off",
spatial = "on",
time = "year",
reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0611529616991593.
mcod_l2_v2 <- run_extra_optimization(mcod_l2, nlminb_loops = 1, newton_loops = 1)sanity(mcod_l2_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large# Large cod model
mcod_l3 <- sdmTMB(density_cod_large ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + breakpt(phi_sc),
data = cod_cpue_large,
mesh = spde_large_cod,
family = tweedie(link = "log"),
spatiotemporal = "off",
spatial = "on",
time = "year",
reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0155203509221983.
mcod_l3_v2 <- run_extra_optimization(mcod_l3, nlminb_loops = 1, newton_loops = 1)sanity(mcod_l3_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large# Large cod model
mcod_l4 <- sdmTMB(density_cod_large ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + s(phi_sc),
data = cod_cpue_large,
mesh = spde_large_cod,
family = tweedie(link = "log"),
spatiotemporal = "off",
spatial = "on",
time = "year",
reml = FALSE)
mcod_l4_v2 <- run_extra_optimization(mcod_l4, nlminb_loops = 1, newton_loops = 1)sanity(mcod_l4_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large#mcmc_res_lc <- residuals(mcod_l_v2, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
#qqnorm(mcmc_res_lc);qqline(mcmc_res_lc)Small flounder
# Small flounder model
# mfle_s <- sdmTMB(density_fle_small ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) +
# breakpt(oxy_sc) + temp_sc,
# data = fle_cpue_small,
# mesh = spde_small_fle,
# family = tweedie(link = "log"),
# spatiotemporal = "off",
# spatial = "on",
# time = "year",
# reml = FALSE)
#
# mfle_s_v2 <- run_extra_optimization(mfle_s, nlminb_loops = 1, newton_loops = 1)# Small flounder model
# mfle_s2 <- sdmTMB(density_fle_small ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) +
# breakpt(oxy_sc) + temp_sc,
# data = fle_cpue_small,
# mesh = spde_small_fle,
# family = tweedie(link = "log"),
# spatiotemporal = "off",
# spatial = "on",
# time = "year",
# reml = FALSE)
#
# mfle_s2_v2 <- run_extra_optimization(mfle_s2, nlminb_loops = 1, newton_loops = 1)#mcmc_res_sf <- residuals(mfle_s_v2, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
#qqnorm(mcmc_res_sf);qqline(mcmc_res_sf)Large flounder
# Large flounder model
# mfle_l <- sdmTMB(density_fle_large ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) +
# breakpt(oxy_sc) + temp_sc,
# data = fle_cpue_large,
# mesh = spde_large_fle,
# family = tweedie(link = "log"),
# spatiotemporal = "off",
# spatial = "on",
# time = "year",
# reml = FALSE)
#
# mfle_l_v2 <- run_extra_optimization(mfle_l, nlminb_loops = 1, newton_loops = 1)# Large flounder model
# mfle_l2 <- sdmTMB(density_fle_large ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) +
# oxy_sc*temp_sc,
# data = fle_cpue_large,
# mesh = spde_large_fle,
# family = tweedie(link = "log"),
# spatiotemporal = "off",
# spatial = "on",
# time = "year",
# reml = FALSE)
#
# mfle_l2_v2 <- run_extra_optimization(mfle_l2, nlminb_loops = 1, newton_loops = 1)#mcmc_res_lf <- residuals(mfle_l_v2, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
#qqnorm(mcmc_res_lf);qqline(mcmc_res_lf)Results
AIC
AIC(mcod_s_v2, mcod_s1_v2, mcod_s2_v2, mcod_s3_v2, mcod_s4_v2)
#> df AIC
#> mcod_s_v2 36 81998.42
#> mcod_s1_v2 37 81908.64
#> mcod_s2_v2 37 81832.75
#> mcod_s3_v2 36 82128.92
#> mcod_s4_v2 36 81727.10
AIC(mcod_l_v2, mcod_l1_v2, mcod_l2_v2, mcod_l3_v2, mcod_l4_v2)
#> df AIC
#> mcod_l_v2 36 110363.0
#> mcod_l1_v2 37 110283.7
#> mcod_l2_v2 37 110213.7
#> mcod_l3_v2 36 110344.2
#> mcod_l4_v2 36 110023.9For both small and large cod, the smooth phi model is best
ggplot(cod_cpue_large, aes(phi, density_cod_large)) + geom_point() + stat_smooth()
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
ggplot(cod_cpue_small, aes(phi, density_cod_small)) + geom_point() + stat_smooth()
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
tidy(mcod_l4_v2)
#> term estimate std.error
#> 1 as.factor(quarter)1 4.42736002 0.4252533
#> 2 as.factor(quarter)4 4.44623473 0.4251376
#> 3 as.factor(year)1994 0.17447088 0.1251899
#> 4 as.factor(year)1995 0.36223149 0.1256474
#> 5 as.factor(year)1996 0.49711011 0.1235093
#> 6 as.factor(year)1997 -0.52187125 0.1276621
#> 7 as.factor(year)1998 -0.64041354 0.1249065
#> 8 as.factor(year)1999 -0.57771772 0.1255800
#> 9 as.factor(year)2000 -0.33946949 0.1319292
#> 10 as.factor(year)2001 -0.25852717 0.1196450
#> 11 as.factor(year)2002 0.46512739 0.1191739
#> 12 as.factor(year)2003 -0.10063961 0.1235872
#> 13 as.factor(year)2004 -0.42438061 0.1231739
#> 14 as.factor(year)2005 0.22425024 0.1165817
#> 15 as.factor(year)2006 0.32941648 0.1173735
#> 16 as.factor(year)2007 0.33068053 0.1143035
#> 17 as.factor(year)2008 0.65553792 0.1138592
#> 18 as.factor(year)2009 0.77452186 0.1116785
#> 19 as.factor(year)2010 1.20435709 0.1121921
#> 20 as.factor(year)2011 0.76709423 0.1130294
#> 21 as.factor(year)2012 0.51482769 0.1158932
#> 22 as.factor(year)2013 0.21750426 0.1141133
#> 23 as.factor(year)2014 0.15737085 0.1158155
#> 24 as.factor(year)2015 0.35017962 0.1137032
#> 25 as.factor(year)2016 0.18424103 0.1133089
#> 26 as.factor(year)2017 -0.05126881 0.1139330
#> 27 as.factor(year)2018 -0.37601593 0.1156271
#> 28 as.factor(year)2019 -0.72588172 0.1387845
tidy(mcod_s4_v2)
#> term estimate std.error
#> 1 as.factor(quarter)1 2.85042323 0.3406691
#> 2 as.factor(quarter)4 2.63170504 0.3412020
#> 3 as.factor(year)1994 -0.69689551 0.1593284
#> 4 as.factor(year)1995 0.12906713 0.1539344
#> 5 as.factor(year)1996 -0.29763762 0.1581533
#> 6 as.factor(year)1997 -0.86644518 0.1553678
#> 7 as.factor(year)1998 -0.44999953 0.1470866
#> 8 as.factor(year)1999 -0.45410484 0.1492065
#> 9 as.factor(year)2000 -0.13818808 0.1546355
#> 10 as.factor(year)2001 0.26872266 0.1408705
#> 11 as.factor(year)2002 0.66548736 0.1423063
#> 12 as.factor(year)2003 -0.49700848 0.1527587
#> 13 as.factor(year)2004 0.29538027 0.1437326
#> 14 as.factor(year)2005 0.69840669 0.1375207
#> 15 as.factor(year)2006 0.05817676 0.1439762
#> 16 as.factor(year)2007 0.37674494 0.1384376
#> 17 as.factor(year)2008 0.54291810 0.1375249
#> 18 as.factor(year)2009 0.54721732 0.1363738
#> 19 as.factor(year)2010 0.79860508 0.1388352
#> 20 as.factor(year)2011 0.36139637 0.1395769
#> 21 as.factor(year)2012 0.28997785 0.1413632
#> 22 as.factor(year)2013 1.04934351 0.1340687
#> 23 as.factor(year)2014 0.77760213 0.1368409
#> 24 as.factor(year)2015 0.43911386 0.1377468
#> 25 as.factor(year)2016 0.24854388 0.1372006
#> 26 as.factor(year)2017 0.27336163 0.1360151
#> 27 as.factor(year)2018 0.02334050 0.1381793
#> 28 as.factor(year)2019 0.22814807 0.1583841
large_pars <- tidy(mcod_l2_v2,
effects = "fixed", conf.int = TRUE) %>%
filter(!grepl('year', term)) %>%
filter(!grepl('quarter', term)) %>%
mutate(size_cl = "Adult")
#> filter: removed 26 rows (84%), 5 rows remaining
#> filter: removed 2 rows (40%), 3 rows remaining
#> mutate: new variable 'size_cl' (character) with one unique value and 0% NA
small_pars <- tidy(mcod_s2_v2,
effects = "fixed", conf.int = TRUE) %>%
filter(!grepl('year', term)) %>%
filter(!grepl('quarter', term)) %>%
mutate(size_cl = "Juvenile")
#> filter: removed 26 rows (84%), 5 rows remaining
#> filter: removed 2 rows (40%), 3 rows remaining
#> mutate: new variable 'size_cl' (character) with one unique value and 0% NA
coefs <- bind_rows(large_pars, small_pars)
coefs <- coefs %>%
mutate(term2 = term,
term2 = ifelse(term == "temp_sc", "Temperature", term2),
term2 = ifelse(term == "oxy_sc", "Oxygen", term2),
term2 = ifelse(term == "depth_sc", "Depth", term2),
term2 = ifelse(term == "oxy_sc:temp_sc", "Temperature x Oxygen", term2)
)
#> mutate: new variable 'term2' (character) with 3 unique values and 0% NA
coefs <- coefs %>%
mutate(term2 = term,
term2 = ifelse(term == "temp_sc", "Temperature", term2),
term2 = ifelse(term == "oxy_sc", "Oxygen", term2),
term2 = ifelse(term == "depth_sc", "Depth", term2),
term2 = ifelse(term == "oxy_sc:temp_sc", "Temperature x Oxygen", term2)
)
#> mutate: no changes
# Coefficients
ggplot(coefs, aes(term2, estimate, color = size_cl)) +
geom_hline(yintercept = 0, linetype = 2, color = "gray50", size = 0.5) +
geom_point(size = 3, position = position_dodge(width = 0.2)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = 1,
position = position_dodge(width = 0.2)) +
labs(x = "Predictor", y = "Standardized coefficient") +
theme(legend.position = "right",
axis.text.x = element_text(angle = 90)) +
NULL 
2. Marginal effects plot
Weird that it’s unimodal and declines at high phi’s. Maybe super coastal?
# Marginal
nd <- data.frame(
phi_sc = seq(min(cod_cpue_small$phi_sc) + 0.5,
max(cod_cpue_small$phi_sc) - 0.2, length.out = 100),
year = 2015L,
quarter = 1,
oxy_sc = mean(cod_cpue_small$oxy_sc),
temp_sc = mean(cod_cpue_small$temp_sc),
depth_sc = mean(cod_cpue_small$depth_sc)
)
p_adu <- predict(mcod_l4_v2, newdata = nd, se_fit = TRUE, re_form = NA, xy_cols = c("X", "Y"))
p_juv <- predict(mcod_s4_v2, newdata = nd, se_fit = TRUE, re_form = NA, xy_cols = c("X", "Y"))
p <- bind_rows(p_adu %>% mutate(size_cl = "Adult"),
p_juv %>% mutate(size_cl = "Juvenile")) %>%
mutate(phi = (phi_sc * sd_phi) + mean_phi)
#> mutate: new variable 'size_cl' (character) with one unique value and 0% NA
#> mutate: new variable 'size_cl' (character) with one unique value and 0% NA
#> mutate: new variable 'phi' (double) with 100 unique values and 0% NA
ggplot(p, aes(phi, exp(est), color = size_cl, fill = size_cl,
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_line() +
geom_ribbon(alpha = 0.4, color = NA) +
geom_vline(xintercept = 1, linetype = 2, color = "grey30") +
coord_cartesian(expand = 0) +
scale_x_continuous(breaks = seq(-2, 40, by = 2))
- Predictions on maps
pred_grid
#> # A tibble: 249,453 × 14
#> X Y depth year oxy temp lon lat sub_div ices_…¹ depth…² T_K
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl>
#> 1 306 6012 13.2 1993 4.11 9.78 12.0 54.2 24 37G2 21 283.
#> 2 306 6012 13.2 1994 6.27 9.66 12.0 54.2 24 37G2 21 283.
#> 3 306 6012 13.2 1995 5.42 9.63 12.0 54.2 24 37G2 21 283.
#> 4 306 6012 13.2 1996 6.11 9.99 12.0 54.2 24 37G2 21 283.
#> 5 306 6012 13.2 1997 5.78 9.84 12.0 54.2 24 37G2 21 283.
#> 6 306 6012 13.2 1998 6.45 8.56 12.0 54.2 24 37G2 21 282.
#> 7 306 6012 13.2 1999 5.53 10.2 12.0 54.2 24 37G2 21 283.
#> 8 306 6012 13.2 2000 4.37 12.5 12.0 54.2 24 37G2 21 286.
#> 9 306 6012 13.2 2001 5.81 9.91 12.0 54.2 24 37G2 21 283.
#> 10 306 6012 13.2 2002 5.26 9.92 12.0 54.2 24 37G2 21 283.
#> # … with 249,443 more rows, 2 more variables: oxy2 <dbl>, phi <dbl>, and
#> # abbreviated variable names ¹​ices_rect, ²​depth_sd
plot_map_fc +
geom_raster(data = pred_grid, aes(X*1000, Y*1000, fill = phi)) +
facet_wrap(~year) +
scale_fill_viridis()
#> Warning: Removed 9882 rows containing missing values (geom_raster).
- weighted phi by size group and year
pred_grid <- pred_grid %>%
mutate(quarter = factor(4),
depth_sc = (depth - mean(depth)) / sd(depth),
phi_sc = (phi - mean(phi)) / sd(phi))
#> mutate: new variable 'quarter' (factor) with one unique value and 0% NA
#> new variable 'depth_sc' (double) with 6,238 unique values and 0% NA
#> new variable 'phi_sc' (double) with 249,453 unique values and 0% NA
large_pred <- predict(mcod_l4_v2, newdata = pred_grid)
small_pred <- predict(mcod_s4_v2, newdata = pred_grid)
# Small cod
wm_phi_small <- small_pred %>%
group_by(year) %>%
summarise("Median" = hutils::weighted_quantile(v = phi, w = exp(est), p = c(0.5)),
"1st quartile" = hutils::weighted_quantile(v = phi, w = exp(est), p = c(0.25)),
"3rd quartile" = hutils::weighted_quantile(v = phi, w = exp(est), p = c(0.75))) %>%
ungroup() %>%
mutate(group = "Juvenile")
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'group' (character) with one unique value and 0% NA
# Large cod
wm_phi_large <- large_pred %>%
group_by(year) %>%
summarise("Median" = hutils::weighted_quantile(v = phi, w = exp(est), p = c(0.5)),
"1st quartile" = hutils::weighted_quantile(v = phi, w = exp(est), p = c(0.25)),
"3rd quartile" = hutils::weighted_quantile(v = phi, w = exp(est), p = c(0.75))) %>%
ungroup() %>%
mutate(group = "Adult")
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'group' (character) with one unique value and 0% NA
env_phi <- pred_grid %>%
group_by(year) %>%
summarise("Median" = quantile(phi, p = c(0.5)),
"1st quartile" = quantile(phi, p = c(0.25)),
"3rd quartile" = quantile(phi, p = c(0.75))) %>%
ungroup() %>%
mutate(group = "Environment")
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'group' (character) with one unique value and 0% NA
phi_series <- bind_rows(wm_phi_small, wm_phi_large, env_phi)
ggplot(phi_series, aes(year, Median, color = group, fill = group)) +
geom_line() +
geom_ribbon(aes(ymin = `1st quartile`, ymax = `3rd quartile`), alpha = 0.2, color = NA) +
facet_wrap(~group) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) +
guides(fill = "none", color = "none") +
labs(y = "phi", x = "Year") +
theme_plot() +
theme(axis.text.x = element_text(angle = 90)) +
NULL
ggplot(phi_series, aes(year, Median, color = group)) +
geom_line() +
scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) +
labs(y = "phi", x = "Year", color = "") +
theme_plot() +
theme(axis.text.x = element_text(angle = 90)) +
NULL